Sampling rare switching events in biochemical networks 
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Bistable biochemical switches are ubiquitous in gene regulatory networks and signal transduction 
pathways. Their switching dynamics, however, are difficult to study directly in experiments or con- 
ventional computer simulations, because switching events are rapid, yet infrequent. We present a 
simulation technique that makes it possible to predict the rate and mechanism of flipping of biochem- 
ical switches. The method uses a series of interfaces in phase space between the two stable steady 
states of the switch to generate transition trajectories in a ratchet-like manner. We demonstrate its 
use by calculating the spontaneous flipping rate of a symmetric model of a genetic switch consisting 
of two mutually repressing genes. The rate constant can be obtained orders of magnitude more 
efficiently than using brute-force simulations. For this model switch, we show that the switching 
mechanism, and consequently the switching rate, depends crucially on whether the binding of one 
regulatory protein to the DNA excludes the binding of the other one. Our technique could also be 
used to study rare events and non- equilibrium processes in soft condensed matter systems. 
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Biochemical switches are essential for the functioning 
of living cells. These switches are networks of chemi- 
cal reactions that exhibit more than one stable steady 
state; in the presence of noise, flipping can occur be- 
tween these states. Well-characterized examples include 
the lysis-lysogeny switch in bacteriophage A Q and the 
lac repressor in E. Coli [^ |j, |^ . Experimental and theo- 
retical studies have established the presence of bistability 
in other biochemical networks, including those regulating 
the cell cycle and developmental fate [M S S H S ll^ ■ 
In addition, synthetic switches have been constructed in 

vivo IHIiilia. 

Computational modeling has an important role to play 
in explaining the properties of biochemical switches. A 
stochastic approach is required to obtain the mechanism 
and rate of switching, since these switches are flipped 
by noise. Examples of such approaches are the chemical 
Langevin technique |l4| , analysis of the chemical master 
equation [15J or stochastic simulation techniques. Sim- 
ulation algorithms that generate trajectories consistent 
with the chemical master equation include the Gillespie 
Algorithm [IE 113 and StochSim [I^. Where spatial res- 
olution is required, methods such as Green's Function 
Reaction Dynamics can be used 19J. 

Biochemical switches are often very difficult or impos- 
sible to simulate using the above techniques in a brute- 
force manner. This is because they can be extremely 
stable, showing few or no flips during the accessible sim- 
ulation time. The average number of spontaneous tran- 
sitions from the lysogenic to the lytic state for bacte- 
riophage A, for example, is about one in 10^ bacterial 
generations |20|, |2l| . New methods are therefore required 
to model such important but rare events in biochemical 
networks. 

Techniques for the simulation of rare events have been 
developed in the field of soft condensed matter physics 



|22| . Recent developments focus on the transition path 
ensemble (TPE). For a rare transition between stable 
states A and B, this is the set of all 'reactive' trajec- 
tories leading from A to B (transition paths). Analysis 
of the TPE gives detailed information on the transition 
mechanism and leads to a prediction of the rate con- 
stant. Transition Path Sampling (TPS) methods have 
been developed to generate members of this ensemble in 
a computationally efficient way [23, |2J| . TPS has been 
applied to a wide variety of problems, including chem- 
ical reactions in solution, conformational transitions in 
biopolymers and protein folding [25j. 

Biochemical switches, however, differ fundamentally 
from these problems. As we shall discuss, in simulations 
of reaction networks the stationary distribution of states 
is generally not known a priori. As a result, TPS meth- 
ods cannot straightforwardly be applied. 

In this article, we present a new scheme for sam- 
pling the TPE and computing the rate constant. This 
"Forward Flux Sampling" (FFS) method is eflicient and 
straightforward. It does not require prior knowledge of 
the phase-space density and can be applied to simulations 
of biochemical networks. The method could also be im- 
plemented in any other stochastic dynamics scheme. To 
our knowledge, FFS constitutes a novel approach to sam- 
pling the TPE. Rather than generating transition paths 
one at a time (as in TPS), a large number of paths are 
grown simultaneously from state A to state _B in a series 
of connected layers. 

As an application of the FFS method, we have calcu- 
lated the spontaneous flipping rate of a simple genetic 
switch, consisting of two mutually repressing genes. We 
show, in agreement with previous work J2y| , that the sta- 
bility of this switch is greatly enhanced when the opera- 
tor regions for the two genes are mutually exclusive, and 
that this is due to an important change in the flipping 



mechanism. 



Rate Expression 



Background 

In this article, the FFS method is used to calculate 
switching rates for biochemical networks simulated with 
the Gillespie Algorithm |16Lll7l| . This algorithm is an ap- 
plication to chemical reactions of the kinetic Monte Carlo 
technique _27J, first introduced by Bortz et al ^8|. The 
system is specified by a set of chemical components {X} 
and a list of allowed reactions, together with their rate 
constants. The concentrations {nx} of all the compo- 
nents are assumed to be homogeneous in space; the state 
of the system at any instant in time is defined by {nx{t)}. 
The concentrations {nx{t)} are propagated stochasti- 
cally in time, assuming each reaction to be a Poisson 
process. This time propagation is consistent with the 
chemical master equation, so that a Gillespie simulation 
is in fact a numerical solution of the master equation. 

An important feature of the Gillespie Algorithm, and 
of other methods for simulating reaction networks, is that 
the distribution of states, i.e. the phase space density, is 
not known a priori, but is an output of the simulation. 
The phase space density can be obtained by solving the 
chemical master equation |29j , but this is generally a de- 
manding task, which is indeed often the motivation for 
carrying out a Gillespie simulation. 

Transition Path Sampling (TPS) has been developed to 
study rare events in condensed-matter systems. In TPS, 
paths belonging to the TPE are obtained by importance 
sampling in trajectory space. New paths connecting sta- 
ble states A and B are generated by making changes to 
existing paths. A new path is accepted or rejected ac- 
cording to its weight in the TPE, which depends on the 
phase space density of its initial point, as well as the 
transition probability for each subsequent step. Without 
prior knowledge of the stationary distribution of states, 
however, this approach cannot conveniently be applied. 

The FFS algorithm which we present in this article 
differs fundamentally from these methods. Rather than 
generating transition paths one at a time, many paths 
are grown simultaneously, in a series of layers, each of 
which forms the basis for the next one. Prior knowledge 
of the stationary distribution of states is not required. 
FFS is well suited for convenient and efficient calcula- 
tion of switching rates in biochemical reaction networks. 
The FFS method is not limited to reaction networks: al- 
though it cannot be used for systems with deterministic 
dynamics, it is applicable to any stochastic dynamical 
scheme, such as Langevin or Brownian Dynamics. In 
this context, it could be used to study rare events in soft 
condensed matter systems such as protein folding and 
crystal nucleation, or non-equilibrium processes such as 
DNA or RNA stretching. 



The expression for the rate constant that is used in the 
FFS algorithm is the same as that described by van Erp et 
al |2d| • The transition occurs between two phase space re- 
gions A and B, which must both be "stable" in the sense 
that if the the system is placed outside these regions, it 
will rapidly evolve in time towards one of them. A and 
B are characterized by the functions Ha^x) and hsix) 
(where x denotes all co-ordinates of the phase space: in 
the case of the Gillespie algorithm the concentrations of 
all the system components), such that: 



hA{x) — liix^A, else hA{x) = 
hsix) = 1 if a; e B, else hsix) = 



(1) 



We also define the functions h^ and h^, which depend 
not only on x{t) but also on the history of the system: 
Ha = 1 if the system was more recently in A than in 
B, and is zero otherwise, while /ig = 1 if the system was 
more recently in B than in A, and is zero otherwise. Thus 
^A + hs = 1 for any point on any path in phase space. 

The rate constant kAB for transitions from region A to 
region B is given by: 



, ^A.B 

KAB = - ' 

Ha 



(2) 



Here, ^a,b is the average number of trajectories per 
unit time entering region B, coming directly from A (i.e. 
which were in A more recently than they were in B). 
Here, the overbar denotes an average over all phase space 
points, with their associated histories. 

The flux ^A,B in Equation |5J is difficult to obtain 
accurately from a simulation because the system makes 
few, if any, spontaneous crossings from A to i? in a typical 
run. To alleviate this problem, a parameter A is chosen, 
such that the functions Ha and hs can be written as: 



hA{x) = 1 if A(a;) < Xa, else hA{x) = 
h-Bix) = 1 if X{x) > Xb, else hsix) = 



(3) 



An increasing series of values of A, {Ai . . . A„}, is then 
chosen, such that Ai > Xa and A„ < Xb. These must 
constitute non-intersecting surfaces in phase space. It is 
not necessary for A to be the reaction coordinate, merely 
that Xa and Xb describe the two stable states (the ex- 
act positioning of these surfaces is not critical). More- 
over, the system should not reach any A^+i before it has 
crossed the preceding surface A^. Defining P{Xi+i\Xi) as 
the probability that a trajectory which passes through 
Ai coming directly from A [i.e. having been in A more 
recently than it last crossed A^), will subsequently reach 
the surface A^+i before returning to A, equation ||2Jl can 
be written as: 



kAi 



$A,_ 



$A, 



hA hA 



-P{Xb\Xi] 



(4) 



Expression Q indicates that the total flux from A to 
B is simply the total flux from A to Xi, multiplied by 
the probability that a trajectory reaching Ai from A will 
eventually arrive in B, before returning to A. A key 
point is that P(Ab|Ai) can be expressed as the product 
of the probabilities of reaching each successive interface 
from the previous one, without returning to A: 



^(AslAi) = H Pi\+i\\) X P(Ab|A, 



(5) 




so that 



A A ^ A,-, A, , A.^ A, o 



FIG. 1: The first stage of the FFS method. 



n-l 



-logP(AB|Ai) = -^logP(A,;+i|A,)-logP(AB|A„) 

(6) 
Expressions ||2J) and Q)-® are used in the FFS method 
to calculate fc^s- 



Forward Flux Sampling 

The first stage of the FFS algorithm involves the choice 
of the parameter A and values for Xa, Xb and {Ai . . . A„}. 
In any one chemical reaction step, the system must be 
able to cross at most one surface Ai. Of course, some 
choices of A will lead to more efficient path sampling than 
others, but we shall demonstrate in the next sections that 
for a typical genetic switch, a rather simple definition of A 
gives very satisfactory results. It is also convenient to de- 
fine a series of "sub-surfaces" {A,- ... X^}, in between 
each pair of surfaces Ai and Ai+i, such that X^ — Xi and 

Aj — Ai+1. 

A simulation is then carried out starting from a point 
in region A. After an equilibration period, the value of 
A is monitored during a run of length T. Whenever the 
trajectory crosses the surface Ai, coming directly from 
A, a counter Nf is incremented. If Nf is less than a 
user-defined number Ci, the phase space co-ordinates of 
the system are also stored. The run is then continued. 
After simulation time T, one is left with a collection of Ci 
points at or just beyond Ai, as well as a measurement of 
the flux ^A,i/hA = Nf/T. This procedure is illustrated 
schematically in Figure Q] crossings of surface Ai that 
are labeled with a black circle contribute to Nf and to 
the collection of points at Ai . |3£,J 

Figure El illustrates the next stage of the algorithm. 
The collection of points at Ai is used to initiate a large 
number Afi of short simulation trial runs. In each of these 
trials, a phase space point from the collection at Ai is 
chosen at random. This is then used as the starting point 
for a simulation run, which is continued until the system 
crosses either A2 or A^. During this run, the maximum 
value of A, Xmax, achieved by the system is recorded. 
Counters N( for all the sub-surfaces A-^"' < Xmax are then 




FIG. 2: The second stage of the FFS method. 



incremented by one. After Mi trials, a good estimate is 
obtained for P{X['^\Xi) = N{/Mi, for 1 < j < mi. We 
note that P(aJ^'|Ai) - 1 and P(Ai'"'^|Ai) = P(A2|Ai). 

During the trial procedure outlined above, one also 
makes a new collection of C2 points at or just beyond 
the surface A2: these are the final phase space points of 
those trial runs starting from Ai which make it to A2 . The 
number M2 of trials must be large enough to generate C2 
points at A2 . The values of M2 , C2 and Mi and Ci for all 
the subsequent surfaces 2 < i < n are chosen by the user: 
the d should be large enough to allow good sampling of 
the phase space. 

The trial run procedure is repeated for each subsequent 
surface Ai , starting from the collection of Ci phase space 
points generated by the successful runs from Ai_i. Even- 
tually Xb is reached, and one is left with a series of his- 
tograms P(A^^]^|Ai), for 1 < i < n and 1 < j < rrii. Us- 
ing equation Q , these histograms can be fitted together 
to obtain a smooth curve P(A|Ai) [Sllll^l, the value of 
which at A = Ab is P(A_b|Ai). The rate constant kAB is 
obtained on multiplying P(A_b|Ai) by the flux ^As/hA 
calculated in the flrst stage of the algorithm. 

It is important to remark that the FFS algorithm does 
not assume that the distribution of phase space points at 
the interfaces {Ai . . . A„} is equal to the stationary dis- 
tribution of states. For the example which we present in 
this paper, this turns out to have significant consequences 
for the transition mechanism. 



Application: A Genetic Switch 

We have applied the FFS method to a simphfied model 
of a genetic toggle switch 26, 33, 34] . This model could 
be regarded as a minimal representation of the lysis- 
lysogeny switch in bacteriophage A [l| ; a synthetic switch 
of this type has also been constructed in vivo [ll| . 

The model switch consists of two proteins A and B and 
their corresponding genes A and B. A and B form ho- 
modimers A2 and B2 which can bind to the DNA strand 
(here labeled O) and influence transcription. When the 
dimer A2 is bound to the DNA, gene B is not tran- 
scribed, while B2, when bound, correspondingly blocks 
transcription of gene A: thus A and B mutually repress 
one another's production. Both proteins are also de- 
graded in the monomer form. We consider two versions 
of this switch: the "general" switch, in which both dimers 
can bind simultaneously to the DNA, forming the species 
OA2B2, and the "exclusive" switch, in which only one 
dimer can be bound at any time. The exclusive switch 
models the case where the operator regions of genes A 
and B are overlapping. 

The switch is represented by the set of reactions H7a(l . 



2A^A2 

O + A2 ^ OA2 

OA2 + B2 ?^ OA2B2 

0^0 + A 

OA2 -^ OA2 + A 

A^0 



2B ^ B2 (7a) 

O + B2 ^ OB2 (7b) 

OB2 + A2 ^ OA2B2 (7c) 

O ^ O + B (7d) 

OB2 -^ OB2 + B (7e) 

B ^ (7f) 



The asterisk indicates that reaction (|7c|l happens only 
for the general switch. Here, we study a symmetrical ver- 
sion of the switch: the rate constants for the reactions on 
the left and right-hand sides of scheme H7a() are identical. 
These are all expressed in terms of the protein produc- 
tion rate constant k (for reactions (|7d|l and ifTeji l. so that 
the unit of time in our calculations is k~^ . The rate con- 
stants for both the forward and backward dimerization 
reactions H7a|) are 5fc. Binding to the DNA occurs with 
rate constant 5fe and dissociation of the complex with 
rate constant k (reactions (|7b|l and lf7c|l l. Finally, the 
rate constant for protein degradation, reactions IjTqI . is 
0.25A;. These parameters are chosen such that in a simu- 
lation using the Gillespie algorithm, the switch flips be- 
tween the A- and B-rich states at a rate that can be mea- 
sured by brute-force simulation. This allows us to test 
the FFS method. The model is, of course, highly sim- 
plistic: our aim is here to demonstrate the FFS scheme 
using a simple example. 




FIG. 3: A as a function of time (in units of fe~^) for a typical 
simulation run, for the general (a) and exclusive (b) switches. 
Insets show probability -P(A) of observing a particular value 
of A, calculated over a total simulation time of 1 x lO^fc"^ 
(general switch, (a)) and 5 x lO^fc"^ (exclusive switch, (b)). 



Results 

Figure |3| shows the results of Gillespie simulations of 
the model genetic switch, in the general (a) and exclusive 
(b) cases. The difference A in the total copy numbers of 
the two proteins, A = Nb — Na, is plotted as a function 
of time, where Na and Nb are defined by: 



Na = ?^A + 2nA2 + 2noA2 
Nb = JT-B + 2nB2 + 2?iob2 



2noA2B2 
2noA2B2 



(8) 



Of course, noA2B2 = for the exclusive switch. Noting 
the different scales on the time axis, the exclusive switch 
(b) has a much lower flipping rate than the general switch 
(a) , in agreement with previous work J2y| . The probabil- 
ity P(A) of obtaining a particular value of A is shown in 
the insets, demonstrating clearly that both the general 
and exclusive switches are bistable. 

Using long brute-force Gillespie simulations, the rate 
constant for the transition from the A-rich to the B-rich 
state was obtained for each switch. We define phase space 
region A to be where A < —25, and region B to be 
where A > 25. The system flips stochastically between 
the state where hj\^ = 1 and hjs = {i.e. it was most 
recently in A) and that where /i^ — and hg = I (it 
was most recently in B). The times t between flipping 
events are distributed according to a Poisson distribu- 
tion p(t) = fcyiB exp [— fc^si], where fc^s = ksA since 
the switch is symmetrical in A and B. The rate con- 
stant kAB can conveniently be measured by fitting the 
cumulative distribution F{t) = f^ dt' p{t') to the func- 
tion 1 — exp [—kAst]. This procedure resulted in values 
of kAB — (4.21 ± 0.05) X 10"^fc for the general switch 
and kAB = (9.4 ± 0.2) x 10~^A; for the exclusive switch 
(using simulation runs of total length 3 x lO^fc"^ [12546 
flips observed] and 9 x lO^fc^^ [8808 flips observed] re- 
spectively). 

We next re-calculated kAB using FFS. The surfaces 
{Ai...A„}, were defined in terms of A: i.e. A = A. 
Regions A and B are given by A = A < A^ and 
A = A > As, respectively. To be sure that the exact 
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FIG. 4: P(A|Ai) for the general (top) and exclusive (bottom) 
switches, where Xa = —24 and As = 25. Black lines: brute- 
force simulation results (averaged over a total time 1 x 10*fc~^ 
for the general switch and 9 x lO^fc"^ for the exclusive switch) ; 
Red lines: FFS results, averaged over 10 independent calcu- 
lations. 



values of A^, As and {Ai . . . A„} did not affect the result, 
these parameters were varied in a series of separate calcu- 
lations. In all cases, we set Ai = A^, and at each surface 
Aj, C, = 10000 points were stored and M, = 100000 
shooting trials were made. Results were averaged over 
10 independent calculations, to obtain error bars similar 
to those of the brute-force results. 

Figure^shows P(A|Ai) as a function of A, for the gen- 
eral and exclusive switches. This function can also be ob- 
tained from an analysis of the brute-force simulation tra- 
jectories: the figure shows excellent agreement between 
the brute- force results (shown in black) , and those of FFS 
(shown in red). 

Tablel^hsts the results for / = $a,i//ix and P(Ab|Ai), 
as well as kAB, for various choices of Xa and As. The 
number n of surfaces used in each calculation is also 
listed. In all cases, the rate constant fc^s is in good 
agreement with the brute-force simulation result. An ad- 
ditional test is provided by measuring / and P(A|Ai) us- 
ing the brute-force runs, to provide a second brute-force 
estimate for kAB- This value is also listed in table |l] and, 
as expected, is in good agreement both with the value ob- 
tained by the fitting to Poisson statistics, and with the 
FFS resuhs. 

Table HTl shows the relative CPU time required for each 
of the calculations. Even for the general switch with a 
relatively fast flipping rate, estimates offers are obtained 
by the FFS method 3 — 6 times faster than using brute- 
force simulations, with similar accuracy. In the case of 
the exclusive switch, the FFS method is 40 — 90 times 
more efficient. Table HTl also demonstrates that the CPU 
time required for an FFS calculation does not increase 
as the rate kAB decreases, in contrast to the time for 
brute-force calculations. Thus FFS should allow the cal- 
culation of rates of even very rare flipping events within 



General switch 






A_g n 


f/k X 10"^ 


P(As|Ai) X 10-3 


kAs/k X 10-5 


30 14 


2.97 ±0.01 


1.41 ±0.03 


4.19 ±0.07 


25 11 


1.33 ±0.01 


3.10 ±0.06 


4.11 ±0.07 


20 9 


0.392 ± 0.003 


10.5 ±0.1 


4.13 ±0.04 


25 - 


1.339 ±0.004 


3.15 ±0.05 


4.22 ±0.06 
4.21 ±0.05 


Exclusive switch 






As n 


f/k X 10-2 


P(As|Ai) X 10-5 


kAs/k X 10-'^ 


30 16 


2.98 ±0.01 


3.2 ±0.1 


9.5 ±0.3 


25 11 


1.211 ±0.007 


7.8 ±0.3 


9.5 ±0.3 


20 10 


0.282 ± 0.002 


33.7 ±0.8 


9.6 ±0.2 


25 - 


1.2112 ±0.0004 


7.70 ±0.09 


9.3 ±0.1 

9.4 ±0.2 



TABLE I: Results for / = $a,i//i^, P(Ab|Ai) and kAB, for 
the general and exclusive switches. FFS results are averaged 
over 10 independent runs, using n surfaces, as described in 
the text. Brute-force results, in bold type, are averaged over 
runs of length 3x lO^fe-^ (general) and 9x lO^fc-^ (exclusive). 
The upper value of kAB is calculated using equation Q a^nd 
the lower value using the fitting of F{t) described in the text. 



reasonable computational time. As an example, we have 
calculated the rate of flipping for a more stable version 
of the exclusive switch, in which the rate constant for 
protein degradation is reduced to 0.175A:. Using the FFS 
method, we obtain a rate kAB — (1-92 ± 0.09) x lO^^fc, 
500 times slower than the switch considered above. This 
result would have been extremely difficult to obtain using 
brute-force simulation. 

The results presented in tableQlshow that for the same 
mean number of protein molecules, the exclusive switch 
has a hipping rate approximately 50 times slower than 
that of the general switch, in agreement with previous 
work |2a|. In order to elucidate the origin of this dif- 
ference, we have analysed the hipping mechanism. The 
FFS method generates a collection of switching trajecto- 



FFS 


Brute-Force 


As 


CPU 


CPU 


general 20 


1.9 


general 6.8 


general 25 


1.1 




general 30 


1.3 




exclusive 20 


2.1 


exclusive 90.2 


exclusive 25 


1 




exclusive 30 


1.7 





TABLE II: Relative CPU time required to calculate the val- 
ues of kAB given in tabled] using parameter values as in the 
text. For the FFS calculations, the calculation of the flux 
/ = '^A,i/hA accounted for 20 - 40% of the total CPU time. 
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FIG. 5: Five randomly chosen transition paths, plotted as 
a function of Na and Nb (a): general switch (b): exclusive 
switch. 



ries (members of the TPE). Figure shows a sample of 
five of these transition paths, for the general switch (a) 
and the exclusive switch (b), plotted as a function of Na 
and Nb. To obtain these paths, we begin with the col- 
lection of partial trajectories that reach Xb from A„ and 
trace these back via the intervening surfaces to Xa- It 
is clear from figure [3 that in the general switch, protein 
A is lost before protein B is gained, so that the transi- 
tion passes through a region of phase space where both 
Na and Nb are low. However, this is not the case for 
the exclusive switch. An important quantity associated 
with transition paths is the committor, Pb{X). This 
is the probability that a new simulation trajectory fired 
from point x will reach region B before A |2^. We have 
measured Pb{x) for points along the trajectories in the 
transition path ensembles generated using FFS, for the 
general and exclusive switches. FigureEfa) shows Pb, as 
well as A = (rtB + 2nB2 ) — (rtA + 2riA2 ) ^^d the occupancy 
of the operator sites, as functions of time for typical tran- 
sition paths. A measures the difference in the number of 
free protein molecules: it is similar but not identical to A. 
A key point is that for the general switch (figure E^(i)), 
the operator makes two important changes of state, from 
OA2 to OA2B2 early in the transition process, and later 
from OA2B2 to OB2. Both of these changes influence 
Pb- For the exclusive switch (figure|Hli(ii)), however, the 
operator is intermittently in states OA2 and OB2 during 
the transition. 

In figure El (b), values of A and operator occupancies 
are shown, averaged over the paths in the TPE, as a 
function of the committor Pb ■ Results are shown for the 
general (i) and exclusive (ii) switches. Figure El (c) anal- 
yses the state points along the paths in the TPE which 
have values of Pb = 0.2, Pb = 0.5 and Pb = 0.8. For 
each of these values of the committor, points are grouped 
according to their operator state. For each Pb and oper- 
ator state, the histograms in figure Efc) show the prob- 
ability distribution function p(A). Clearly, for points on 
a constant Pb surface, the operator state and the num- 
ber of free molecules are correlated: when A is bound to 
the DNA, on average, a larger A is required to obtain 
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FIG. 6: (a): Typical transition paths, for the general (i) and 
exclusive (ii) switches. A = (hb + 2nB2) — ("-a + 2nA2), 
Pb and operator occupancy (color-coded) are plotted versus 
time. The path for the general switch (a) happens to be 
longer than that for the exclusive switch (b), but this is not 
true of all paths in the ensemble, (b): Average values of A 
and average operator occupancies, as functions of the com- 
mittor Pb, for the general (i) and exclusive (ii) switches, (c): 
Probability distribution p(A), for selected values of Pb- p(A) 
is divided into color-coded contributions from the different 
operator states. The sum of the areas of the histograms for 
a particular Pb is unity. For all plots, 100 test trajectories 
were used to estimate Pb- 



a particular value of Pb than when B is bound. This 
shows that Pg, and hence the reaction co-ordinate, de- 
pends not only upon the difference in copy number A, 
but also on the state of the operator. Importantly, the 
plots of figure (b) and (c) are not symmetric on mak- 
ing the transformations Pb — > 1 — Pb, A -^ —A and 
OA2 ^^ OB2. This demonstrates that the distribution of 
transition paths does not follow the steady state phase 
space density, which is symmetric on interchangin g A and 
B, since the switch is by construction symmetric |35j- It 
also means that the TPE for the reverse transition, from 
B to A, would occupy a different region of phase space, 
as compared to the TPE for the transition from A to 
B (see also supplementary material) . The origin of this 
asymmetry is that the dynamics of our system involves ir- 



reversible reactions (Reactions (J7dl7i^ ). Indeed, our sys- 
tem does not satisfy microscopic reversibility. Finally, 
figure El clearly demonstrates why the elimination of the 
operator state OA2B2 for the exclusive switches enhances 
its stability with respect to that of the general switch. In 
the general switch, as soon as a B2 dimer is produced by 
some rare fluctuation, it can bind to the DNA, switch off 
the production of A and thereby accelerate the flipping 
of the switch. For the exclusive switch, however, any B2 
dimer that is produced must wait for a second fluctua- 
tion by which A2 is released from the DNA, before it can 
bind. This is the origin of the enhanced stability of the 
exclusive switch. 



to understand the origin of the enhanced stability of the 
exclusive switch. The FFS method will be easily appli- 
cable to many important biochemical switches, for which 
prediction of the rates and pathways of switching should 
lead to a better understanding of the design principles 
underlying their stability. 

We thank Peter Bolhuis for very helpful discussions 
and Daan Frenkel, Rutger Hermsen and Harald Tepper 
for their careful reading of the manuscript. The work is 
part of the research program of the "Stichting voor Fun- 
damenteel Onderzoek der Materie (FOM)", which is fi- 
nancially supported by the "Nederlandse organisatie voor 
Wetenschappelijk Onderzoek (NWO)". 



Discussion 



This article presents the Forward Flux Sampling (FFS) 
method for the calculation of the rates of rare events in 
stochastic kinetic simulation schemes such as those used 
for biochemical reaction networks. In contrast to pre- 
viously developed methods for sampling the transition 
path ensemble [2J,|33, FFS does not require knowledge 
of the phase space density. It is this feature that makes 
it possible to study rare events in biochemical networks. 
Our algorithm samples the TPE in a way that is, to our 
knowledge, new: many paths are grown simultaneously 
from state A to state B in a series of layers of partial 
paths, each layer forming the basis for the next. The 
phase space separating stable states A and B is traversed 
by the algorithm in a "ratchet-like" manner, making the 
method highly suitable for very rare events, where one- 
at-a-time path generation tends to be inefficient. FFS is 
not applicable to systems whose dynamics is determinis- 
tic. It could be used, however, in combination with any 
stochastic simulation technique. This will make it useful 
for a wide range of problems in soft condensed matter 
systems, including rare events and non-equilibrium pro- 
cesses. 

We have demonstrated our method using stochastic 
simulations of a simple genetic switch consisting of two 
mutually repressing genes. Following earlier work |2fil |. 
we compare the case where both protein products can 
bind simultaneously as dimers to the DNA (the general 
switch), to that where each protein dimer excludes the 
binding of the other (the exclusive switch). The results 
obtained using FFS are in good agreement with those 
of long brute- force simulations for both switches. The 
computational time required for the FFS calculations is 
far less than for the brute-force simulations, and in addi- 
tion, does not increase as the rate constant decreases. In- 
deed, using FFS we could simulate a switch that was too 
stable to be studied using brute-force calculations. By 
analysing the transition path ensembles we were able to 
discover the differences between the flipping mechanisms 
for the general and exclusive switches. These allow us 
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Supplementary information: asymmetric transition paths 
for a symmetric genetic switch 
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Figure 1: Average operator occupancy as a function of A = (ne + SnBa) ~ 
(riA + SnAs), fo^ transition paths of the symmetric toggle switch consisting 
of two genes that mutually repress each other (Eqs. 7 in main text), (a): 
general switch; (b): exclusive switch. The transition paths from A to S, 
as indicated by the solid lines, are asymmetric on replacing ^4 by i? {i.e. 
A -^ —A and OA2 OB2). This is illustrated by the dotted lines, which, 
in fact, correspond to the transition paths from B to A. This shows that the 
transition trajectories from A to i? do not coincide with those from B to A. 
The origin of the asymmetry is that the switch does not obey microscopic 
reversibility. 



